#include "HY1C_out.h"
/* --------------------------------------------------------------- */
/* get_par.c - computes photosynthetically active radiation        */
/*                                                                 */
/* Inputs:                                                         */
/*     l2rec - level-2 structure containing one complete scan      */
/*             after atmospheric correction.                       */
/* Outputs:                                                        */
/*     par   - Photosynthetically Active Radiation just above the  */
/*             surface from SeaWiFS level 1b instantaneous         */
/*             radiances at 412, 443, 490, 510, 555, and 670 nm.   */
/*                                                                 */
/* Algorithm Provided By: Robert Frouin,                           */
/*                        Scripps Institution of Oceanography      */
/*                                                                 */
/* Written By: B. A. Franz, SAIC GSC, SIMBIOS, 30 July 1999        */
/*                                                                 */
/* --------------------------------------------------------------- */

#include <stdlib.h>
#include <math.h>
#include "l12_proto.h"
#include "smi_climatology.h"

int calc_par_seawifs_(
               int32_t   *year, 
               int32_t   *month, 
               int32_t   *day, 
               float  *time, 
               float  *lat, 
               float  *lon, 
      	       float   Rad[], 
               float  *SolZen, 
               float  *ViewZen, 
               float  *RelAz, 
               float  *TauA865, 
               float  *AngstromI, 
               float  *Dobson, 
               float  *wv,
      	       float  *SurfPress, 
               float  *wave, 
               float  *F0, 
               float  *kO3, 
               float  *TauMol, 
               float  *PAR );

int calc_par_modis_( 
               int32_t   *year, 
               int32_t   *month, 
               int32_t   *day, 
               float  *time, 
               float  *lat, 
               float  *lon, 
      	       float   Rad[], 
               float  *SolZen, 
               float  *ViewZen, 
               float  *RelAz, 
               float  *TauA865, 
               float  *AngstromI, 
               float  *Dobson, 
               float  *wv,
      	       float  *SurfPress, 
               float  *wave, 
               float  *F0, 
               float  *kO3, 
               float  *TauMol, 
               float  *PAR );


int calc_par_viirs_( 
               int32_t   *year, 
               int32_t   *month, 
               int32_t   *day, 
               float  *time, 
               float  *lat, 
               float  *lon, 
      	       float   Rad[], 
               float  *SolZen, 
               float  *ViewZen, 
               float  *RelAz, 
               float  *TauA865, 
               float  *AngstromI, 
               float  *Dobson, 
               float  *wv,
      	       float  *SurfPress, 
               float  *wave, 
               float  *F0, 
               float  *kO3, 
               float  *TauMol, 
               float  *PAR );


void get_par(l2str *l2rec, float par[])
{
    int32_t  ip, ib, ipb, iw;

    float Lt[NBANDS];

    float angst;
    float taua;
    float delphi;
    float ozone;

    int32_t  year  = *l2rec->year;
    int32_t  day   = *l2rec->day;
    float hour  = *l2rec->msec/3.6E6;

    static int32_t lastDay  = -999;
    static int32_t month    = -999;
    static int32_t monthDay = -999;

    static int32_t mask = SEAICE | LAND | HIGLINT;

    static float lambda[NBANDS];     
    static float Fobar [NBANDS];     
    static float Taur  [NBANDS];     
    static float kO3   [NBANDS];  
    static int   bindx [NBANDS];   
    static int   nbands;
    static int   firstCall = TRUE;

    if (firstCall) {

        firstCall = FALSE;

        /* Initialize climatologies */
        smi_climatology_init(l2rec->input->alphafile,day,ALPHA510);
        smi_climatology_init(l2rec->input->tauafile ,day,TAUA865 );

        
        switch (l2rec->sensorID) {

	  case HMODISA:
          case HMODIST:
            nbands = 3;
            bindx[0]=2; bindx[1]=6; bindx[2]=7;
            break; 
	  case SEAWIFS:
            nbands = 6; 
            for (ib=0; ib<nbands; ib++)
                bindx[ib]=ib;
            break;
	  case VIIRS:
            nbands = 5; 
            for (ib=0; ib<nbands; ib++)
                bindx[ib]=ib;
            break;
          default:
            HY1C_out("PAR not supported for this sensor.\n");
            exit(1);
            break;
	}

        /* Get band-pass dependent quantities */
        for (iw=0; iw<nbands; iw++) {
            ib = bindx[iw];
            lambda[iw] = l2rec->fwave[ib];
            kO3   [iw] = l2rec->k_oz[ib]/1000.0;
            Fobar [iw] = l2rec->Fobar[ib];
            Taur  [iw] = l2rec->Tau_r[ib];
        }
    }

    /* PAR function wants month and day-of-month */
    if (*l2rec->day != lastDay) {
        int16 mon;
        int16 mday;
        yd2md(year,day,&mon,&mday);
        month    = mon;
        monthDay = mday;
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        /* Grab radiances for this pixel */
        for (ib=0; ib<nbands; ib++) {
            ipb = ip*NBANDS + bindx[ib];
            Lt[ib] = l2rec->Lt[ipb];
        }

        /* Skip pixel if masked */
        if (Lt[0] <= 0.0 || (l2rec->flags[ip] & mask) != 0 || l2rec->solz[ip] > 83.0) {
	  //HY1C_out("%f %d %f\n",Lt[0],l2rec->flags[ip],l2rec->solz[ip]);
            par[ip] = BAD_FLT;
            l2rec->flags[ip] |= PRODFAIL;
            continue;
        }        

        /* Want relative azimuth in 0 - 180.0 */
        delphi = l2rec->sena[ip] - l2rec->sola[ip];
        if (delphi >  180.0) delphi -= 360.0;
        if (delphi < -180.0) delphi += 360.0;

        /* Want ozone in Dobson units */
        ozone = l2rec->oz[ip]*1000;

        /* Get angstrom and AOT from climatology */
        angst = smi_climatology(l2rec->lon[ip],l2rec->lat[ip],ALPHA510);
        taua  = smi_climatology(l2rec->lon[ip],l2rec->lat[ip],TAUA865);

        switch (l2rec->sensorID) {

	  case HMODISA:
          case HMODIST:

//            calc_par_modis_(
//                  &year,
//                  &month,
//                  &monthDay,
//                  &hour,
//                  &l2rec->lat[ip],
//                  &l2rec->lon[ip],
//                  Lt,
//                  &l2rec->solz[ip],
//                  &l2rec->senz[ip],
//                  &delphi,
//                  &taua,
//                  &angst,
//                  &ozone,
//                  &l2rec->wv[ip],
//                  &l2rec->pr[ip],
//                  lambda,
//                  Fobar,
//                  kO3,
//                  Taur,
//                  &par[ip]);
	    break;

          case SEAWIFS:

//            calc_par_seawifs_(
//                  &year,
//                  &month,
//                  &monthDay,
//                  &hour,
//                  &l2rec->lat[ip],
//                  &l2rec->lon[ip],
//                  Lt,
//                  &l2rec->solz[ip],
//                  &l2rec->senz[ip],
//                  &delphi,
//                  &taua,
//                  &angst,
//                  &ozone,
//                  &l2rec->wv[ip],
//                  &l2rec->pr[ip],
//                  lambda,
//                  Fobar,
//                  kO3,
//                  Taur,
//                  &par[ip]);
            break;

          case VIIRS:

//            calc_par_viirs_(
//                  &year,
//                  &month,
//                  &monthDay,
//                  &hour,
//                  &l2rec->lat[ip],
//                  &l2rec->lon[ip],
//                  Lt,
//                  &l2rec->solz[ip],
//                  &l2rec->senz[ip],
//                  &delphi,
//                  &taua,
//                  &angst,
//                  &ozone,
//                  &l2rec->wv[ip],
//                  &l2rec->pr[ip],
//                  lambda,
//                  Fobar,
//                  kO3,
//                  Taur,
//                  &par[ip]);
            break;

          default:
            HY1C_out("PAR not supported for this sensor.\n");
            exit(1);
            break;
	}

        /* Convert to E/D/m^2 */
        par[ip] *= 1.193;

    }
}



/*
Subject: 
       PAR routine
   Date: 
       Fri, 23 Jul 1999 08:55:29 -0700
  From: 
       Robert Frouin <rfrouin@ucsd.edu>
    To: 
       chuck@seawifs, gfargion@simbios, gene@seawifs, wang@simbios, franz@seawifs
    CC: 
       jmcpherson@ucsd.edu




Greetings:

A routine to compute daily PAR from SeaWiFS level 1b radiances is available
at the following address: http://genius.ucsd.edu/~john/SeaWiFS_dir/ under
the rubrique "PAR subroutine and test program".

The routine requires as input year, month, day, time, latitude, longitude,
SeaWiFS radiances in the first 6 spectral bands, solar zenith angle,
viewing zenith angle, relative azimuth angle, aerosol optical thickness at
865 nm, Angstrom coefficient, ozone amount, and surface pressure. Routine
output is daily PAR.

Thus a daily PAR value is computed for each instantaneous SeaWiFS
observation, clear or cloudy. Diurnal variations are taken into account
statistically. The algorithm is described succintly in the routine.

During our discussion at GSFC in June, a first routine was supposed to be
developed to provide a normalized cloud/surface albedo and then a second
routine to compute daily PAR from the normalized albedo, and the second
routine was to be applied when binning to the 9 km resolution. Now daily
PAR is obtained using a single routine, which is more convenient.

The binning to the 9 km resolution should be done as follows. First,
weight-average the daily PAR estimates obtained from all SeaWiFS
observations during the same day at each location (there might be several
SeaWiFS observations of a surface target during the same day). The weight
is the cosine of the sun zenith angle for the SeaWiFS observation. That is:

        PAR_avg = sum{cos[tetas(i)]*PAR(i)}/sum{cos[tetas(i)]}

Second, simply average the values at all the locations within the 9 km bins.

The routine requires aerosol data, ozone amount, surface pressure. If these
parameters are missing (-999 or less), default values are used. If Eric
Vermote's aerosol climatology (or Menghua's) is not available yet, please
use default values for tests.

At this time, the statistical diurnal function does not depend on latitude,
longitude, and date, but will depend on these parameters in the second
version of the code. Creating a date and location dependent function
requires analysing several years of ISCCP data. We have the data, but a
couple of weeks is needed to accomplish the task.

I will present the algorithm at the SeaWiFS atmospheric correction meeting
next week, and prepare a detailed documentation.

Best, Robert.


Robert Frouin
Scripps Institution of Oceanography
University of California San Diego
9500 Gilman Drive
La Jolla, CA 92093-0221
Voice Tel.: 619/534-6243
Fax Tel.: 619/534-7452


*/
